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Abstract 

Applying the mathematical circulation theory of Markov chains, we investigate 
the synchronized stochastic dynamics of a discrete network model of yeast cell-cycle 
regulation where stochasticity has been kept rather than being averaged out. By com- 
paring the network dynamics of the stochastic model with its corresponding determin- 
istic network counterpart, we show that the synchronized dynamics can be soundly 
characterized by a dominant circulation in the stochastic model, which is the natural 
generalization of the deterministic limit cycle in the deterministic system. Moreover, 
the period of the main peak in the power spectrum, which is in common use to char- 
acterize the synchronized dynamics, perfectly corresponds to the number of states in 
the main cycle with dominant circulation. Such a large separation in the magnitude 
of the circulations, between a dominant, main cycle and the rest, gives rise to the 
stochastic synchronization phenomenon. 
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I Introduction 

Synchronization is an important characteristics of many biological networks [32\ [35] whose 
dynamics has been modelled traditionally by deterministic, coupled nonlinear ordinary dif- 
ferential equations in terms of regulatory mechanisms and kinetic parameters [2C1, 4j . Two 
important classes of biological networks which have attracted wide attentions in recent 
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years are neural networks [3T] and cellular biochemical networks [8]. A deterministic 
model, however, only describes the averaged behavior of a system based on large popula- 
tions; it can not capture the temporal fluctuations of a small biological system with either 
extrinsic or intrinsic noise. For example, neuronal firing has an inherent variability, and 
many biochemical regulatory networks inside a cell consist of molecular species with low 
concentrations, where stochastic models with chemical master equations (CME) based on 
biochemical reaction stoichiometry, molecular numbers, and kinetic rate constants should 
be developed. Such an approach has already provided important insights and quantitative 
characterizations of a wide range of biochemical systems |19t [33], (H [2"4"1 [57] . We refer to 
[34j for a good introduction to stochastic modelling in biology. 

As there is a growing awareness and interest in studying the effects of noise in bi- 
ological networks, it becomes more and more important to quantitatively characterize 
the synchronized dynamics mathematically in stochastic models, because the concepts of 
limit cycle and fixed phase difference no longer holds in this case. Instead, physicists and 
biologists always have to characterize synchronized dynamics by the distinct peak of its 
power spectrum or just only by observing the stochastic trajectories, which however may 
cause ambiguities in the conclusion. Therefore, a logical generalization of limit cycle in 
stochastic models needs to be developed. 

On the other hand, from the view of statistical physics, these stochastic models for 
systems cell biology exhibit nonequilibrium steady states (NESS) in which nonequilibrium 
cycle fluxes necessarily emerge [22]. The relation between an NESS and traditional nonlin- 
ear dynamics, in particular bistability [25] and limit-cycle oscillations, has been discussed 
recently [231124]. 

We have recently developed a rather complete mathematical theory of nonequilibrium 
steady state (NESS, p3j) of open systems, taking Markov chains as the model of biochem- 
ical systems. One of the most important concepts in the mathematical theory of NESS 
is the circulation(also called cycle fluxes) [H], which corresponds to the cycle kinetics in 
open chemical systems |10[ 122], Actually, since the stochastic circulation in NESS is de- 
fined as the time-averaged frequency of each cycle in the sense of trajectory, it can also be 
regarded as a quantitative characterization for the intensities of stochastic synchronized 
dynamics in different attractor basins. This paper will thoroughly investigate the interplay 
described above in a currently interested stochastic cell-cycle network model. 

Cell cycle is one of important biological processes whose underlying molecular networks 
have been extensively studied. Tyson and his coworkers have studied mathematical models 
for yeast cell cycle based on differential equations [2]. Their key finding is that cell cycle 
is a hysteresis loop, via saddle-node bifurcations, driven by the periodic changes in cell 
mass due to cell growth and division [4j. The biological "check points" correspond to the 



steady states of the dynamical system [3J. 

However, in many cellular biochemical modelling, a detailed, molecular number-based 
CME model is not warranted because of a lack of quantitative experimental data. Thus 
alternatively, one may try to develop discrete state network, such as Hopfield network and 
Boltzmann machine, model of a complex biological system using the available information 
on the activation and repression from one signaling molecules to another, e.g., the kind of 
signaling wiring diagram like Fig. [TJ 




Figure 1: The cell-cycle network of the budding yeast. Each node represents a protein or 
a protein complex. Arrows are positive regulation, T-lines are negative regulation, dotted 
T-loops are degradation. It is just Fig. 1 from [38] with permission. 

Recently, Li, et.al. [18] have developed a discrete deterministic Boolean model by 
applying the approach of Hopfield for neural networks to the yeast cell-cycle regulatory 
network with 11 nodes according to Fig. [TJ They have found that the topology of the 
network provides significant robustness of the dynamics toward the check points, i.e., 
steady states. The deterministic Boolean model has been further extended to incorporate 
stochastic dynamics [38J. Several biologically interesting results have been obtained; these 
include the numerical evidence for the existence of a single dominant cell cycle and its 
robustness under a large range of noise level. 

In the present paper, we will investigate the synchronized stochastic dynamics of the 
discrete network model of yeast cell-cycle regulation where stochasticity has been kept 
rather than being averaged out from the start, applying the mathematical circulation 
theory of Markov chains. By comparing the network dynamics of the stochastic model 
with its corresponding deterministic network counterpart, we show that the synchronized 
dynamics can be soundly characterized by a dominant circulation in the stochastic model, 



which is the natural generalization of the deterministic limit cycle in the deterministic 
system. 

It is shown that a large separation in the magnitude of the circulations, between a dom- 
inant, main cycle and the rest, gives rise to the stochastic synchronization phenomenon, 
and the power spectrum of the trajectory has a main peak, whose period converges per- 
fectly to the number of states in the dominant cycle. Furthermore, the net circulation 
of the dominant cycle increases monotonically with the noise-strength parameter 0, ap- 
proaching to its deterministic limit. Together, these observations provide a clear picture 
of the nature of the synchronization in a stochastic network in terms of the probabilistic 
circulation of NESS. The main part of the present paper, Section IV presents our study on 
the synchronization and dominant circulation in the cell-cycle network through numerical 
computations. 

For the completeness of the work, we first give a theoretical sketch of some relevant 
results on biological networks in Section II, including a classification of the deterministic 
and stochastic Boolean networks and their correspondence. Then a short review of the 
mathematical theory of stochastic circulation for Markov chains is introduced and applied 
to Hopfield network and Boltzmann machine in Section III. It is shown that the stochastic 
Boolean network is reversible if and only if the matrix T in the model is symmetric, and 
the net NESS circulation is strictly positive as long as the probability of the directed cycle 
is larger than that of its reversal. Finally, the conclusion and some further discussions are 
provided in Section V. 

II Theoretical Sketch of Some Relevant Results on Biolog- 
ical Networks: Hopfield networks, Boltzmann machines, 
and Their Relations 

We first give a brief account of the deterministic Hopfield networks and its stochastic in- 
carnation, the Boltzmann machines. Since the influential work of J.J. Hopfield in 1980s' 

[T2] , the deterministic Boolean (Hopfield) network has been applied to various fields 
of sciences. Amit [1] has introduced a temperature-like parameter (3 that characterizes 
the noise in the network and constructed a probabilistic Boolean network called Boltz- 
mann machine. The Hopfield networks and Bolzmann machines have found wide range of 
applications in biology. They can be mainly categorized into several classes as follows pQ. 

Let us suppose TV is a fixed integer, S = {1, 2, • • • , iV}. We take the state space as 
{-M} s - 

Model A: deterministic 

Al (discrete time, synchronous, McCulloch-Pitts) : Denote the state of the nth step as 



X n = (Xn, X\, • • • , X~"), then the dynamic is 

at \ N 



Xi +1 = signj [ ]T TijXi -Ui\,ttJ2 TijXl ± U t ; (1) 



and if ^2j = iTijXn = Ui, then X % n+1 randomly choose 1 or —1 with probability | respec- 
tively. Ui (1 < i < iV), given a priori, are called the threshold of the zth unit. 

A2 (continuous time, synchronous): Every state has an exponentially distributed 
stochastic waiting-time, with mean waiting-time A" 1 , then chooses the next state by the 
same rule of model A\. 

A3 (discrete time, asynchronous, Hopfield): The neurons are updated one by one, 
in some prescribed sequence, or in a random order. If the previous state a satisfies 
YljLi TijPj > Ui, then <Tj changes to be 1, otherwise changes to be —1. 

A4 (continuous time, asynchronous): Every state has an exponentially distributed 
stochastic waiting-time, with rate constant A, then chooses the next state by the same 
rule of model A3. 

Remark II. 1 The (—1, 1) representation and the (0, 1) representation are equivalent. Bi- 
ologists usually use the (0, 1) representation in modelling biological networks, where and 
1 denote the resting and activated states, respectively. Physicists usually use the (—1, 1) 
representation for spin systems. 

Remark II. 2 Note that by deterministic, we mean the transition from one state to the 
next is deterministic. But the systems with continuous-time will still behave stochastically 
due to the Poisson nature in the transition time. 

Model B: stochastic Boltzmann machines 

Bl (discrete time, synchronous): Consider the Markov chain 

{X n = {X l n ,X 2 n ,--- ,X*),n = 0,1,2,- ■■} (2) 

on state space {— 1, with transition probability given as follows: for each pair of states 
o~,n E { — 1, 1} S , the probability transiting from a to n 

= t~t eMPr]i(Ilf=i T ijVj ~ Ui) (3) 

PaV l\ exp(/3(Ef = i -Ui))+ exp(-/3(Ef =1 T l3 a 3 -Utf' 

where (5{> 0), Ty and Ui (1 < i,j < N) are parameters of the model. 



1 The function sign(a;) = < 



1 x > 0; 
x = 0; 
-1 a;<0. 



B2 (continuous time, synchronous): Consider the continuous-time Markov chain {£ t : 
t > 0} on state space {—1, 1} S . Every state waits an exponential time with meantime A -1 
until choosing the next state by the rule of model Bl. So the transition density matrix is 

q* v = ^P* v , Va,JjG{-l,l} S . (4) 

B3 (discrete time, asynchronous): Denote a % to be the new state which changes the 
sign of the ith coordinate of a. Consider the Markov chain {X n = (X^, • • • , X^), n = 
0, 1, 2, • • • } on state space {—1, 1} S ■ The neurons are updated one by one, in some pre- 
scribed sequence, or in a random order. Then choose the next state according to the 
probability: 



exp(-/^(Ef =1 T^ _£/.)) ( 
~ exp(/?(£f=i TijO-j - Ui)) + exp(-/?(£f =1 T %J a 3 - U)) ' ° " ; ' " 1 



where (3 > 0; and p arj = 0, if n ^ a 1 for each i. 

B4 (continuous time, asynchronous): Every state waits an exponential time with mean- 
time A -1 until choosing the next state by the rule of model 53. So the transition density 
matrix is 

W = *JW> Va € (5) 

The third class given below is a variant of the model Bl. It is included here since it is 
the model used for the probabilistic Boolean network of cell-cycle regulation. 
Model C: deformation 

Fix a > 0. Consider a new Markov chain {X n = (X^, X%, • • • , X^),n = 0, 1, 2, • • • } 
on the state space { — 1,1} taking the model Bl as defined initially, with transition 
probability given as follows: 

N 

P(X n+1 \X n ) = l[P i (Xi +1 \X n ), (6) 



i=l 



where we define 



Pi{X l n+1 \X n ) — 



/ j j = 1 ?>3 n / % 



cxp(/3(Ef=i T lJ ^-(7 l ))+cxp(- / 3(E J / Li TijXl-Ui)) ^r- 

if Ef=i TijXi = Ui and X* n+1 = X* 



l+e- a 

e- a :v V- V 



if V =1 T iJ X n = Ui and X* n+1 = l-Xl 



Remark II. 3 The model C differs from Bl when Ylj=i TijXh — Ui = 0, and the latter is 
also a special case of the former when a = 0. 



Ill Mathematical Circulation Theory of Network Nonequi- 
librium Steady States 



There are many different approaches to the theory of nonequilibrium statistical mechanics 
in the past [91 1 2 1[ ITS], mathematical theories of which have emerged in the last two decades, 
and Jiang et.al. have summarized their results of this theory in a recent monograph [14]. 
The most important concepts in the theory are (i) reversibility of a stationary process that 
corresponds to thermodynamic equilibrium, and (ii) the circulation in a stationary process 
which corresponds to NESS. A key result of the theory is the circulation decomposition of 



III.l Circulation theory of nonequilibrium steady states 

Hill [10] constructed a theoretical framework for discussions of vivid metabolic systems, 
such as active transport, muscle contractions, etc. The basic method of his framework is 
diagram calculation for the cycle flux on the metabolic cycles of those systems [10] . He 
successively found that the result from diagram calculation agrees with the data of the 
numbers of completing different cycles given by random test (Monte Carlo test), but did 
not yet prove that the former is just the circulation rate in the sense of trajectory of a 
corresponding Markov chain. In Chapter 1 and 2 of [14] . Markov chains with discrete 
time and continuous time parameter are used as models of Hill's theory on circulation in 
biochemical systems. The circulation rate is defined in the sense of trajectories and the 
expressions of circulation rate are calculated which coincide with Hill's result obtained from 
diagrams. Hence the authors verify that Hill's cycle flux is equivalent to the circulation 
rate defined in the sense of trajectories. 

Below we only state the circulation theory of discrete-time Markov chains, and refer to 
Chapter 2] for the quite similar circulation theory of continuous-time Markov chains. 

First of all, we state the main results, rewritten in terms of the cycle representation of 
stationary homogeneous Markov chains ([17]. Theorem 1.3.1), which is analogous to the 
Kirchoff's current law and circuit theory in the networks of master equation systems [30] , 

Theorem III.l Given a finite oriented graph G = (V,E) and the weight on every edge 
{uj e > : e € E}. If the weight satisfies the balance equation 



then there exists a positive function defined on the oriented cycles {uj c : c = (ei, e2, • • • , e&), k G 
Z + }, such that 



NESS. 




(7) 




(8) 



For a finite stationary Markov chain X with finite state space S = {1,2,- •• ,N}, 
transition matrix P = {pij : i,j£ S} and invariant distribution ff = {tt\, 7T2, ■ • • , ttn }, one 
can take its state space to be the group of vertexes, and E = {e : e + = i, e~ = j,Pij > 0} 
with weight {u> e = TTiPij}. Then notice that ([7]) is satisfied, since YljPij = 1 f° r each i and 

X! ^iPy = = ^2 KjPjU v * G V, 

we can conclude that 

Corollary III. 2 (cycle decomposition) For an arbitrary finite stationary Markov chain, 

there exists a positive function defined on the group of oriented circuits {oj c : c = 12, ■ ■ ■ , ik), k G 

Z + } such that 

nPij = ^u c J c (i,j), Vi,j G V, (9) 

c 

where J c (i,j) is defined to be 1 if the cycle c includes the path i — > j, otherwise 0. 
Definition III. 3 u; c zs called the circulation along cycle c. 
For any i,j e S,i^ j, 

KiPij ~ TTjPji = ^(wc - w c _)J c (i,j), (10) 
c 

where c_ denotes the reversed cycle of c. Equation (I10p is called the circulation decom- 
position of the stationary Markov chain X. 

It can be proved that generally the circulation decomposition is not unique, i.e. it is 
possible to find another set of cycles C and weights on these cycles {w c \c E C} which also 

fit crop. 

However, the most reasonable choice of circulation definition is the one defined in the 
sense of trajectories form the probabilistic point of view. Along almost every sample path, 
the Markov chain generates an infinite sequence of cycles, and if we discard every cycle 
when it is completed and at the meantime record it down, then we can count the number 
of times that a specific cycle c is formed by time t, which we denote by w Ct t(u>). 

The following theorem is recapitulated from [141 Theorem 1.3.3]. 

Theorem III. 4 Let C n (u>), n = 0, 1,2, ••• , be the class of all cycles occurring until n 
along the sample path {X[(u>)}. Then the sequence (C n iuS) , w c^nitd) / n) of sample weighted 
cycles associated with the chain X converges almost surely to a class (Coo,u; c ), that is, 

Coo = hm C n (cA a.e. (11) 

n— >+oo 

w c = lim c,n ^ - , a.e. (12) 
rw+oo n 



Furthermore, for any directed cycle c = ■ ■ ■ ,i s ) € Coo; the weight w c is given by 

D({ii,i2,--- ,is} c ) i „v 

w c — VixnVi-iiz * * * Pi s -iisPisii ^ D({j} c ) ' 

where D = {dij} = I — P = {5ij — pij} and D(H) denotes the determinant of D with 

0, % + 



rows and columns indexed in the index set H. The function 5ij 
known Kronecker delta function. 



is the well 
1, i = 3, 



It is important to emphasize that the circulations defined in the above theorem also 
satisfy the circulation decomposition relation (llOp , 

The above theorem not only rigorously confirms the Hill's theory, but also gives a prior 
substitute method of the widely used diagrammatic method. The complexity of directed 
diagrams and cycles increases rapidly with the number of states in the model, while the 
determinant interpretation is much more systematic and easy to be applied using the 
mathematics softwares. But it will still cost excessive time to compute these determinants 
if there are hundreds of states in the model. 

Luckily, a Monte Carlo method using the so-called derived chain method [14, Sec- 
tion 1.2] to compute the cycle fluxes has already been developed according to the above 
theorem, the main idea of which is just simply to discard every cycle when it is completed 
and at the meantime record it down so as to count the number of times that a specific 
cycle c is formed by time t(i.e. w c j(u))). Then when the time t is long enough, one gets 
the approximated circulation of the cycle c (i.e. w c » Wc 'l^ ). 

Actually, the Boolean yeast cell-cycle network model discussed in the next section has 
2048 states and we find that the Monte Carlo method is much more efficient than the 
method of determinant interpretation, because the latter is even impossible to be applied 
to such a large model upon a normal computer. 

The relationship between circulation and NESS is as follows, which is recapitulated 
from |14| Theorem 1.4.8]. 

Theorem III. 5 Suppose that X is an irreducible and positive-recurrent stationary Markov 
chain with the countable state space S, the transition matrix P = (pij)ij^s an d the in- 
variant probability distribution U = (n^i^s, Q-nd let {w c : c E Coo} be the circulation 
distribution of X, then the following statements are equivalent: 

(i) The Markov chain X is reversible. 

(ii) The Markov chain X is in detailed balance, that is, 

KiPij = iTjPji,Vi,j € S. 
(Hi) The transition probability of X satisfies the Kolmogorov cyclic condition: 



Vi\iiViiiz ' ' ' Pi a -iisPish ~ PhisPisis-i ' ' ' PhizP^hi 



for any directed cycle c = ■ ■ ■ ,i s ). 

(iv) The components of the circulation distribution of X satisfy the symmetry condition: 



Wc = Wc- , Vc € C, 



'OO- 



Consequently, when the system is in a nonequilibrium steady state, there exists at least 
one cycle, containing at least three states, round which the circulation rates of one direction 
and its opposite direction are asymmetric (unequal), so as to cause a net circulation on 
the cycle. In theoretic analysis, if there is a large separation in the magnitude of the 
circulation, between few dominant, main cycles and the rest, it gives rise to the stochastic 
synchronization phenomenon and helps to distinguish the most important main biological 
pathways, which can be observed in experiments. 

III. 2 Applied to Hopfield network and Boltzmann machine 

The keys to understand synchronization behavior in stochastic models are (i) establishing 
a correspondence between a stochastic dynamics and its deterministic counterpart; and 
(ii) identifying the cyclic motion in the stochastic models. 

In the framework of the stochastic theory, deterministic models are simply the limits 
of stochastic processes with vanishing noise. This is best illustrated in the following 
proposition. 

Proposition III. 6 With the same initial distribution, when (5 — > oo, the model con- 
verges to the model in distribution, for k = 1,2,3,4. 

Proof: From 



We shall further discuss the necessary condition for stochastic Boolean networks to 
have cyclic motion. In the theory of neural networks, one of Hopfield's key results 112] 
is that for symmetric matrix Tij, the network has an energy function. In the theory of 
Markov processes, having a potential function (Kolmogorov cyclic condition in Theorem 
IIII.5|) is a sufficient and necessary condition for a reversible process. A connection is 
established in the following theorem. 

Theorem III. 7 For k = 1,2,3,4, the Markov chain {X n } in the model B^ is reversible 
if and only if Tij = Tji, Vz, j = 1, 2, • • • ,N, i.e. the matrix T is symmetric. 



P->oo eMPiELi Tijaj - Ui)) + exp(-/3(V /V =1 T^aj - Ui)) 




(14) 



Proof: According to theorem IIII.51 we can use the Kolmogorov cyclic condition to com- 
plete the proof. Since every two states are communicative, one only needs to consider the 
cycles consist of three different states. 

Denote these states by a = (a%, a.2, ■ • ' , &n), o = (ci, o~2, • ■ • , &n)j an d V = (^l) V2, • • • > Vn)- 
Consider the cycle a — > a — > n — > a then 



PaaPcrrjPrja 

ParjPricrPcra 
N 



^(P<Ti(Y$Ll T i3 a j - U i)) 



n 



i exp(/3(V = i Tijoij - Ui)) + exp(-/3(X\- =1 l^oy - U { )) 



N 

n 



exp((3r)iCE,f=i TijUj - Ui)) 



i exp(/3(E -Li T^i - Ui)) + eM-P(E7=i T H°i ~ U i)) 



N 

n 



exp(/3a i (Y l f =1 Ti j r]j - Ui)) 



i exp(/3(^^Li Tijr/j - Ui)) + exp(-/3(X)jLi T^j - Ui)) 



' N 

n 



exp(/3r/ i (^;^ 1 Ijj-aj - £/*)) 



i exp(/3(^^Li I^-ay - C/i)) + exp(-P(Y$Li T^-ay - E/ f )) 



exp(/3a i (X;^L 1 lyJfe - 



n 



i exp(/3(X;^Li 3ij»7j - Ci)) + exp(-P(T,j=i TijVj ~ Ui)) 



N 



N 

n 



exp(Pai(52f=i Tijdj - Ui)) 



i=1 exp(/3(X) j= i 2<iOj- - Ui)) + exp(-/3(^ i= i 2ij0-j - Ui)) _ 
exp < Tij[(otiOj + OiT}j + r/ia-,) - (aycr* + <7j77j + njai)] 



exp 



icrj + o-jr/j + ruaj) 



(15) 



necessity: Let rjk = —1, VA; = 1, 2, ■ • • , TV, = —1, V7c ^ i, on = 1, and a k = — 1, VA; 7^ j, 
°j = !• So Ylk,i( T ki - Tik){a k ai + a k f]i + m a i) = 2 (Tij ~ T ji), if the Markov chain 
is reversible then by Theorem IIII.5I T iii) PaaPa-qP-qa = ParfPrjaPaa, we can conclude that 
T ij = T ji , V*,j = l,2,.-- ,JV. 

sufficiency: If T^- = Tji, Vi, j = 1, 2, • • • JV, then from (fT5|) one has p a(T Par]Pria = Pa-qPr/aPaa, Va, a, 77. 



From the above two conclusions, it is obvious that 

Corollary III. 8 When Tij = Tji,Vi,j = 1,2, ••• ,N, there doesn't exist any limit cycle 
consist of more than two states in model A k , k = 1, 2, 3, 4. 



Remark III. 9 Part of the proof for the corollary \III.8\ can be found in the literature, 
such as [36]. The proof we give here, however, provides a new point of view from that of 
stochastic convergence, and is complete. 



Lemma III. 10 For k = 1,2,3,4, the Markov chain {X n } in the model is reversible 
when (3 is zero, and the Markov chain {X n } in the model C is reversible when (5 and a 
are both zero. 

Proof: In the case when j3 and a are both zero, then for arbitrary two states a and 77, 
we have p ar) = |. So the Markov chain {X n } in the model B^, k = 1, 2, 3,4, and model C 
satisfies the Kolmogorov cyclic condition. | 

IV Synchronized Dynamics in a Stochastic Network of Yeast 
Cell-Cycle Regulation 

From biochemical perspective, the microscopic variables for a cellular regulatory network 
are the concentrations, or numbers, of various mRNAs, regulatory proteins, and cofac- 
tors. If all the biochemistry were known, then the dynamics of such a network would be 
represented by a chemical master equation [19] . Unfortunately, much of the required 
information is not available, nor such a "fully-detailed" model will always be useful. Phe- 
nomenologically the concentrations of key players of a biochemical regulatory network 
can often be reduced to two or three states, such as resting state, activated state, inacti- 
vated state, etc. The interactions between these states are usually determined from 
experimental data. 

IV. 1 The stochastic network model of Zhang et. al. 

We now turn to the stochastic model of cell-cycle network as developed by Li et.al. [18J 
and Zhang et.al. |38[ I39j . The model in |18j is a deformation of model Ay with the 
(0, 1) representation instead of the (—1, 1) representation: Denote the state of n-th step 
as X n = (X*, X%, ■ • • , X^), then the dynamic is 

(N \ 1 N 

J2 T iJ X i) + 1 > if'^'/'-,-V/, ••'(): (16) 

while X l n+1 = X l n if J2f=i TijXl = 0. 

The main results of |18| are that the network is both dynamically and structurally 
stable. The biological steady state, known as the Gl phase of a cell cycle, is a global 
attractor of the dynamics; the biological pathway, i.e., the returning to Gl phase after 
perturbation, is in a globally attracting basin. There is no limit cycle in this model. 

Zhang et.al. [38J implemented a probabilistic Boolean network for the cell-cycle protein 
interaction network. They found that both the biological steady state and the biological 
pathway are well preserved under a wide range of noise level. The model in [38] is model 



C, with the trivial difference of taking the (0, 1) representation rather than (—1, 1) repre- 
sentation. 

Notice that, according to Proposition IIII.61 when /3, a — > oo, this model recovers the 
deterministic model in [18] : hence, it is implicit that there is no dominant cycle with 
significant circulation when (3 and a are large. 

There is an important difference between the previous models and the most recently 
developed one in [39J, which will be investigated currently in more detail. The earlier works 
all treat the cell mass (volume), implicitly through protein Cln3 (node 1), as a parameter 
of their models rather than a dynamic part of the models. Under this setting, the Gl 
phase of the cell cycle, (0,0,0,0,1,0,0,0,1,0,0) in our binary representation, is a global 
attractor. The biological cell cycle, however, has a clear sense of cycling with direction. It 
is a time-irreversible process accompanied with positive entropy production; it is a system 
in NESS. The cell mass has been incorporated into the transition dynamics in [39], again 
implicitly characterized by the 7 parameter: When the system is in Gl phase, a signal to 
exit the Gl phase by activating Cln3 (node 1), known as START in cell biology, will come 
with a probability given in Eq. (|17p below. It represents the trigger implicitly associated 
with cell mass. This is the idea of [39], but the method of attack used there is ad hoc with 
no complete picture. 

In the present paper, we focus on the model in [39], thus, is a deformation of model C: 
fix a > 0, 7 > 0, consider the Markov chain {X n = (X^,X^, • • • , X^ ) : n = 0, 1, 2, • • • } 
on the state space {0, 1} S , whose transition probability is as follows: 

N 

P{X n+ i\X n ) = Y\_Pi{X % n+1 \X n ), 



i=l 



where 

Pi{X l n+1 \X n ) 



exp(/3(2X; +1 -l)(Ef =i r^)) 
exp(/3(£f=i TijXi)) + ex P (-/3(£f =1 TyXl)) ' 



if£f =1 ^X^0. 



X l „ 11 — X, 



l-|_ e -o> 1 x — A ni 



]_j_ e - Q ' n+1 ^ v n' 

Xt.,1 = 1 — XL, 



if E =1 T H x n = and X n ^ (0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0) or i > 2. 



Pi(.Y #(+1 |A\,)-< 1+ f Xl (17) 

l_|_ e 7 ' ^n+l 1 n> 

if Ef=i T 1:j Xl = and X n = (0, 0, 0, 0, 1,0, 0, 0, 1, 0, 0). 

The work [39] only introduced this model which is more reasonable than that in [38J, 
but the results and conclusions in [39] are very short and shallow. So the aim of the present 



paper is to give more profound and comprehensive insight into their work, applying the 
mathematical theory of nonequilibrium steady states. 

The matrix T in all the models above [18} [38| [39] according to FigQ] is 
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(18) 



where Ty = 1 for a positive regulation of protein i to protein j and T« = — 1 for a negative 
regulation of protein i to protein j. If the protein i has a self-degradation loop, Ta = —0.1. 

This matrix is not symmetric. Hence, according to theorem IIII.71 the stochastic dy- 
namics has a NESS with nonzero circulation. 

Notice that, when /3, a and 7 — > 00, this model converges to a deterministic model 
similar to that in [18]: Denote the state of n-th step as X n = (X^, X^, ■■■ , X^), then the 
dynamic is 



while X l n+1 = X l n if Ylf=i TijXl = and X n + (0,0,0,0,1,0,0,0,1,0,0) or i > 2; 
Xl +1 = 1 - X\ if Zf=i TijXl = and X n = (0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0). 



The corresponding deterministic model has a limit cycle consist of 13 state. So there 
is also a main cycle whose circulation is dominant when /3, a and 7 are large, which is 
described by Tab. [TJ 

IV. 2 Synchronization and dominant circulation in cell-cycle network 

While Zhang et.al. [381 13T?] have introduced the probabilistic Boolean network model of 
yeast cell cycle and tried to describe its nonequilibrium dynamics, they are not fully aware 
of the important interplay between the the nonequilibrium circulations of this model and 
the observed synchronization phenomenon. 

Furthermore, from numerical computation of the maximal probability flux on cycles, 
the previous work [381 ES] seems to suggest that the cycle with positive circulation is 




(19) 



Cln3 MBF SBF Clnl,2 Cdhl Swi5 Cdc20/Cdcl4 Clb5,6 Sicl Clbl,2 Mcml/SFF Phase 

(1 1 1 0) Start 

(0 11010 010 0) Gi 

(0 11110 010 0) Gi 

(0 1 1 1 0) Gi 

(0 11100 100 0) s 

(0 11100 101 1) G 2 

(0 00100 1 101 1) M 

(0 1 1 1 1) M 

(0 00001 1 Oil 1) M 

(0 00001 1 010 1) M 

(0 1 1 1 1 0) M 

(0 00011 010 0) Gi 

(0 00010 010 0) Gi 



Table 1: The synchronous sequence of 13 states as recorded in Li et al. [18] 

unique. This is not the case. Equipped with Corollary IIII.2I on cycle decomposition , we 
shall show that there are many other cycles with positive circulation in the model, although 
their probability weights are indeed quite small in contrast to the main cycle. This large 
separation in the magnitudes of the weights gives rise to the stochastic synchronization. 
Therefore, this is indeed a characterization of the stochastic synchronized dynamics. 

Indeed, from the theory of nonequilibrium circulation of discrete-time homogeneous 
finite Markov chains introduced in Section III, the net circulation is strictly positive as 
long as the probability of directed cycle is larger than that of its reversal. This is true for 
any cycle on which the Kolmogorov cyclic condition breaks down [14, Theorem 1.4.8], e.g., 
in the model in [39], let states /i = (0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0), a = (1,0, 0, 1, 0, 0, 0, 0, 0, 1,0), 
and v = (0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1). Consider the cycle fi — ► a — *■ v — * /i. Then 

PnaPavPufi _ e 3o+10.9/3 ^ 
PfivPvcrPtrfj, 

Then the net circulation of this cycle is not zero. 

Numerical computations with the famous Gillespie's method [7] are carried out, and 
the results are given in the following figures. The network with 11 binary nodes has a total 
of 2048 number of states. Following Amit (pQ, pp. 75-79), we present the dynamics of the 
network in terms of integers 0, 1, 2, • • • , 2047, where the state (sx, S2, • • ■ , sn) corresponds 
one-to-one to Yll=i s i2 J_1 . This 1-d system obtained is reversible if and only if the 11-d 
system is reversible. 

Fig. [2] is the basic behavior of a random trajectory. The upper panel shows that there 
arises the phenomenon of local rapid synchronization like that observed in [13] during a 
very short time period, when (5 is sufficiently large. The lower panel is a random trajectory 
over a much longer time. 

To further characterize the synchronized dynamics, Fig. [3] shows the Fourier power 
spectrum of the trajectories in Fig. [2] Using MATLAB, the discrete Fourier transform 



Local rapid synchronization 
2000 1 1 — i 1 1 

1500 - 




Figure 2: Stochastic trajectory and synchronization. Simulations are carried out with the 
parameters a = 5, (3 = 6, and 7 = 5. 



for time series {x\, X2, • • • , x n } is defined as 



n 



i(27r/n)(m-l)(fc-l) 
fc=l 



(20) 



which is just the magnitude of frequency — ^-27r, 1 < m < n. 

Therefore, by the Herglotz theorem ([26], P- 331), the power spectrum of discrete 
trajectories has a symmetry y m = y n +2-m- For different sets of parameters, all the calcu- 
lations give the distinguishable main peak in the Fig. [3j It is important to mention that 
different trajectories tend to the same Fig. [3] by ergodicity. 

The single dominant peak implies there exists a global synchronization phenomenon. 
Besides the main peak, the power spectra also have several other smaller peaks, which 
correspond to several cyclic motions with lower magnitudes. Note that the synchronized 
behavior is preserved in representation that maps, one-to-one, from the 11 binary nodes 
to the integers — 2047. It is possible that the map will cause some distortion in the 
power spectrum. To further illustrate the synchronized behavior, Fig. 0] shows the power 
spectra of all the 11 individual nodes in the network. While subtle details are different, 
all exhibit the dominant peak, similar to that of the overall dynamics. This demonstrates 
further the synchronized dynamics in the network. 

Fig. [5] plots the magnitude and the period of the dominant peak of the power spectrum 
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Figure 3: Power spectrum of the overall trajectory with the parameters a = 5, 7 = 5 fixed, 
and different (3: blue: (3=2 A; green: /?=4.8; red: (3 = 6. The discrete Fourier transform 
causes an alias; hence the spectrum is symmetric with respect to ir on the [0, 2ir] interval. 

in Fig. [3J as functions of the noise strength, i.e., the parameter (3. It shows, as we have 
predicted, that the period converges to 13 which corresponds perfectly to the number of 
states in the main cycle when (3 is large. We also put error bars on the upper panel of 
Fig. [5] with various values of (3. 

Finally, Fig. [6] shows how the net circulation of the dominant, main cycle varies with 
(3. It is clearly seen that the net circulation of the main cycle increases monotonically, 
which implies the appearing of more and more distinct synchronization with increasing (3. 
The direction of the net circulation does not change. 

The net circulations of all the cycles are very small when (3, a and 7 are near zero, 
since the system is close to equilibrium state (reversible) according to lemma UlL 101 

And for large (3, the net circulations of all the cycles except the main cycle are very 
small, whose order are about 10 -6 by numerical simulation using the Monte Carlo method 
for circulation of Markov chain according to Theorem IIII.41 All the circulations of non- 
main cycles actually decrease with increasing (3 when (3 is large. 

As in [38], we also notice that there exists an inflection point in the curve in Fig. 
[H This implies a cooperative transition of the net circulation of the main cycle while 
varying the noise level (3. We are currently studying the nature of this phase-transition 
like behavior. The implication of this observation remains to be further elucidated. 

It is also important to point out that the above results are very robust with respect 
to the values of the matrix T and the threshold U, which is similar to the main result in 
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Figure 4: Power spectra of individual nodes show a synchronization among all the nodes 
with the parameters a = 5, 7 = 5 fixed, and different (3: blue: (3=2A; green: /3=4.8; red: 
/?=6. 

The appearance of the peak in the power spectrum is a signature of a coherent reso- 
nance; the existence of circulation means some components of the system are synchronized. 
A mathematical equivalence between the power spectral peaking and the circulation was 
proved by Qian et.al. [29] in the case of Markov chains. It has been generalized to all 
Markov processes by Jiang and Zhang |15j . 

V Conclusion and Discussion 

Quantitative understanding of biological systems and functions from their interacting com- 
ponents presents a significant challenge as well as an unique opportunity for scientists of 
diverse disciplines. Based on the mathematical theory of nonequilibrium steady states, 
especially the cycle circulation distribution, synchronized dynamics in a stochastic yeast 
cell-cycle network is investigated in detail and quantitatively, which is more definite and 
convincing than the previous method often used by physicists and biologists. 
Synchronization and circulation in stochastic network. 

As we know, the occurrence of a deterministic limit cycle in an ordinary differential 
equation(ODE) model is the hallmark of a synchronization phenomenon. For a stochas- 
tic system, such a definite concept no longer holds and a logical generalization needs to 




Figure 5: Magnitude and period of the dominant power-spectral peak as functions of f3, 
with a = 5 and 7 = 5. In the upper panel, the magnitude of the dominant power-spectral 
peak is averaged over 10 simulations. The solid curve is the mean, and the dotted curves 
are the mean ± standard deviation. In the lower panel, the abscissa for the period of the 
cycle, T, is in logarithmic scale. The period approaches to 13 (the horizonal line) with 
increasing (5. 

be investigated. In fact, the very meaning of synchronization requires clarification. In 
this case, the concept of stochastic limit cycle is useful in describing the synchroniza- 
tion phenomenon. For example, in our present model of yeast cell cycle, the trajectory 
concentrates around a main cycle, which can be defined as a stochastic limit cycle, is in 
many aspects very similar to the limit cycle of a deterministic model. Although in many 
situations this concept does not have a clear definition, it has been rigorously discussed 
in the context of Markov chain under the heading of circulation theory for many years 
[27} [28] . We recommend the Chapters 1 and 2 of [14] for the systematic treatment of this 
theory. 

It is shown in the present paper that the stochastic synchronized dynamics in the 
cell-cycle model, first observed in [38], can be nicely characterized in the perspective of 
the theory of nonequilibrium steady states and Markov chain circulations. The stochastic 
Markov chain is reversible (equilibrium state) if and only if the net circulation of every 
cycle vanishes, which is one of the most important theorems in the mathematical theory 
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Figure 6: Net circulation of the main cycle as a function of the noise strength, (3. 
of NESS. 

The stochastic Boolean network is reversible if and only if the matrix T is symmetric, 
and the net NESS circulation emerges as long as the probability of the directed cycle is 
larger than that of its reversal, which is just the preconditions for synchronized dynamics 
to occur. So circulation can be regarded as the quantitative characterization of stochastic 
synchronized dynamics. Those cycles with high circulation perform stronger synchronized 
dynamics, while other cycles with low circulation perform weaker synchronized dynamics. 
Moreover, the distinct global synchronized dynamics can be obviously observed if there is 
a large scale separation in the magnitude of the circulation between an unique main cycle 
and the others. 

On the other hand, circulation as a good characterization of synchronized dynamics 
in stochastic models, corresponds perfectly to the power-spectrum algorithm in determin- 
ing the synchronized dynamics. The magnitudes of the main cycle's circulation and the 
dominant power-spectral peak both increases monotonically with increasing /?, showing 
more and more distinct synchronization in the biological systems. Also the period of the 
dominant power-spectral peak converges to the period(number of states) of the main cycle 
with dominant circulations. 

Stable attractor and robustness in the presence of biological stochasticity. 

Biological systems have to be robust to function in complex (and noisy) environments. 
More robustness could also mean being more evolvable, and thus more likely to survive. 
In the model discussed in the present paper, the unique limit cycle in the deterministic 
network model and the monotonically increasing circulation in the stochastic network 
model both demonstrate the stability of the main stochastic limit cycle. It can be obviously 




seen that the circulation of the main cycle dominate under a large range of noise level 
(Fig. [6|), which consequently means that this biological function is very robust against 
small perturbations. It is directly responsible for this cellular process. 

Moreover, in some stochastic Boolean network model, one can contain both stable 
and unstable stochastic limit cycles, where the unstable stochastic limit cycle surprisingly 
causes vanishing circulation when the noise intensity tends to zero, although it actually 
is a true limit cycle of the corresponding deterministic Boolean network model 
which attracts more than half the states in the state space! This phenomenon, 
called "noise induced global attractive behavior" , is recently observed by us in a stochastic 
Boolean network model of the Siah-l/beta-catenin/pl4/19 ARF loop of the protein p53 
pathways (unpublished work), in which the biological trajectory, starting from the initial 
state where only p53 protein is activated, is attracted to the unique stable stochastic 
limit cycle of the protein states rather than the other unstable one. This fact further 
confirms the robustness of this biological system. Certainly, more implication of our 
observations remain to be further elucidated in the future, when applied to other significant 
bio-problems. 
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